Serratia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
path_serratia <- "../../../../output/2_Oligotyping/2A/Serratia/2A_oligotyping_Serratia_sequences-c1-s1-a0.0-A0-M10"

# load the matrix count table
matrix_count <- read.table(paste0(path_serratia, "/MATRIX-COUNT.txt"), header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP41 S104 S25 S27 S30 S32 S36 S37 S38 S39 S44 S48 S49 S51 S77
A 10 24 10 7 442 3 10 3 4 8 0 26 2 2 0 4 3 2 1200 965 55 2132 46 666 17 2 43 13 314
G 4 6 1 0 207 0 7 2 2 5 1 5 2 1 40 0 5 0 655 298 36 874 7 306 6 0 23 5 0

Taxonomy

# load the fasta table
fasta <- readDNAStringSet(paste0(path_serratia, "/OLIGO-REPRESENTATIVES.fasta"))

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
A A
G G

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## load the reference table
ref_oligo_med2 <- read.table("../../../../output/2_Oligotyping/2B/2B_REF_info_serratia.tsv", sep="\t", header = TRUE)

## select only the 2 oligotypes of Serratia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP41 S104 S25 S27 S30 S32 S36 S37 S38 S39 S44 S48 S49 S51 S77
oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013 10 24 10 7 442 3 10 3 4 8 0 26 2 2 0 4 3 2 1200 965 55 2132 46 666 17 2 43 13 314
oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498 4 6 1 0 207 0 7 2 2 5 1 5 2 1 40 0 5 0 655 298 36 874 7 306 6 0 23 5 0
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013 oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013
oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498 oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498

Metadata

metadata <- read.csv("../../../../metadata/metadata.csv", sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 29 samples ]
## sample_data() Sample Data:       [ 29 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 8511
# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 28 samples ]
## sample_data() Sample Data:       [ 28 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_serratia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Serratia
ps.serratia <- ps.filter %>% subset_taxa(Genus=="Serratia")

# add read depth of Serratia
sample_data(ps.filter)$Read_serratia <- sample_sums(ps.serratia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"        
## [13] "Read_serratia"
sample_data(ps.serratia) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"
# add percent of Serratia
sample_data(ps.filter)$Percent_serratia <- sample_data(ps.filter)$Read_serratia / sample_data(ps.filter)$Read_depth

# round the percent of Serratia at 2 decimals
sample_data(ps.filter)$Percent_serratia <- sample_data(ps.filter)$Percent_serratia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_serratia, Percent_serratia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Strain Field Country Organ Species Run Control Species_italic Strain_italic Read_depth is.neg Read_serratia Percent_serratia
CTC1 CTC1 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 29779 FALSE 14 0.00
CTC10 CTC10 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2609 FALSE 30 0.01
CTC11 CTC11 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 13874 FALSE 11 0.00
CTC12 CTC12 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1146 FALSE 7 0.01
CTC13 CTC13 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 18035 FALSE 649 0.04
CTC14 CTC14 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1708 FALSE 3 0.00
CTC15 CTC15 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 23180 FALSE 17 0.00
CTC2 CTC2 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 30692 FALSE 5 0.00
CTC3 CTC3 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 39920 FALSE 6 0.00
CTC4 CTC4 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2139 FALSE 13 0.01
CTC5 CTC5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 15789 FALSE 1 0.00
CTC6 CTC6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 19753 FALSE 31 0.00
CTC9 CTC9 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 4980 FALSE 3 0.00
NP14 NP14 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 7973 FALSE 0 0.00
NP2 NP2 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 648335 FALSE 0 0.00
NP20 NP20 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 136 FALSE 0 0.00
NP27 NP27 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 1234 FALSE 0 0.00
NP29 NP29 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 203 FALSE 0 0.00
NP30 NP30 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 228 FALSE 0 0.00
NP34 NP34 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 95 FALSE 0 0.00
NP35 NP35 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 196532 FALSE 0 0.00
NP36 NP36 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 249 FALSE 0 0.00
NP37 NP37 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 419340 FALSE 0 0.00
NP38 NP38 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 282479 FALSE 0 0.00
NP39 NP39 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 218684 FALSE 0 0.00
NP41 NP41 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 247152 FALSE 40 0.00
NP42 NP42 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 185157 FALSE 0 0.00
NP43 NP43 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 239335 FALSE 0 0.00
NP44 NP44 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 156879 FALSE 0 0.00
NP5 NP5 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 736159 FALSE 0 0.00
NP8 NP8 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 334799 FALSE 0 0.00
S100 S100 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52486 FALSE 0 0.00
S102 S102 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 3456 FALSE 0 0.00
S104 S104 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52403 FALSE 4 0.00
S105 S105 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55577 FALSE 0 0.00
S106 S106 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 33053 FALSE 0 0.00
S107 S107 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52154 FALSE 0 0.00
S108 S108 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55735 FALSE 0 0.00
S109 S109 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59023 FALSE 0 0.00
S110 S110 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 57377 FALSE 0 0.00
S121 S121 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20361 FALSE 0 0.00
S122 S122 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 9803 FALSE 0 0.00
S123 S123 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20130 FALSE 0 0.00
S124 S124 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 18146 FALSE 0 0.00
S126 S126 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 15235 FALSE 0 0.00
S127 S127 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 24696 FALSE 0 0.00
S128 S128 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 16305 FALSE 0 0.00
S146 S146 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25012 FALSE 0 0.00
S147 S147 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25171 FALSE 0 0.00
S148 S148 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 14164 FALSE 0 0.00
S150 S150 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15081 FALSE 0 0.00
S151 S151 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22944 FALSE 0 0.00
S152 S152 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15082 FALSE 0 0.00
S153 S153 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17040 FALSE 0 0.00
S154 S154 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 9626 FALSE 0 0.00
S160 S160 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 72508 FALSE 0 0.00
S162 S162 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25180 FALSE 0 0.00
S163 S163 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 12333 FALSE 0 0.00
S164 S164 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22368 FALSE 0 0.00
S165 S165 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17731 FALSE 0 0.00
S166 S166 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 13979 FALSE 0 0.00
S167 S167 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 14048 FALSE 0 0.00
S169 S169 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 11553 FALSE 0 0.00
S170 S170 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 8852 FALSE 0 0.00
S18 S18 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4290 FALSE 0 0.00
S19 S19 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44527 FALSE 0 0.00
S20 S20 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42864 FALSE 0 0.00
S21 S21 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33798 FALSE 0 0.00
S22 S22 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 19044 FALSE 0 0.00
S23 S23 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 38172 FALSE 0 0.00
S24 S24 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42355 FALSE 0 0.00
S25 S25 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 47688 FALSE 8 0.00
S26 S26 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 5394 FALSE 0 0.00
S27 S27 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 24558 FALSE 2 0.00
S28 S28 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4503 FALSE 0 0.00
S30 S30 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 25353 FALSE 1855 0.07
S31 S31 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 20417 FALSE 0 0.00
S32 S32 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 12441 FALSE 1263 0.10
S33 S33 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33867 FALSE 0 0.00
S34 S34 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 9367 FALSE 0 0.00
S35 S35 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 11663 FALSE 0 0.00
S36 S36 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33020 FALSE 91 0.00
S37 S37 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 18340 FALSE 3006 0.16
S38 S38 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 54790 FALSE 53 0.00
S39 S39 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 36273 FALSE 972 0.03
S40 S40 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44448 FALSE 0 0.00
S42 S42 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 4107 FALSE 0 0.00
S43 S43 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 9279 FALSE 0 0.00
S44 S44 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 8026 FALSE 23 0.00
S45 S45 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 18150 FALSE 0 0.00
S47 S47 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 1951 FALSE 0 0.00
S48 S48 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56738 FALSE 2 0.00
S49 S49 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 33498 FALSE 66 0.00
S50 S50 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 28481 FALSE 0 0.00
S51 S51 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 61788 FALSE 18 0.00
S52 S52 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 21553 FALSE 0 0.00
S55 S55 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 50447 FALSE 0 0.00
S56 S56 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 42609 FALSE 0 0.00
S57 S57 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49157 FALSE 0 0.00
S58 S58 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 30357 FALSE 0 0.00
S59 S59 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 32798 FALSE 0 0.00
S60 S60 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 44485 FALSE 0 0.00
S61 S61 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49545 FALSE 0 0.00
S63 S63 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53444 FALSE 0 0.00
S64 S64 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 47628 FALSE 0 0.00
S79 S79 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59755 FALSE 0 0.00
S80 S80 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52788 FALSE 0 0.00
S83 S83 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 42272 FALSE 0 0.00
S84 S84 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56676 FALSE 0 0.00
S85 S85 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 41690 FALSE 0 0.00
S86 S86 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 61984 FALSE 0 0.00
S87 S87 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 65958 FALSE 0 0.00
S88 S88 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53102 FALSE 0 0.00
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_serratia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_serratia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=2, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 2 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013.",
##  "oligotype_G (23) | size:2498 / node_N1116 (23) | size:2498.",
data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)
new_names <- c("oligotype_A (27) | size:6013 / N1114 (27) | size:6013.",
               "oligotype_G (23) | size:2498 / N1116 (23) | size:2498."
)

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col_add <- brewer.pal(8, "Accent")

col <- c("oligotype_A (27) | size:6013 / N1114 (27) | size:6013."="#B136F5",
         "oligotype_G (23) | size:2498 / N1116 (23) | size:2498."="#B863FF")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1 taxa and 1 samples ]
## sample_data() Sample Data:       [ 1 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 1 taxa by 1 taxonomic ranks ]
# data pour plot
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_A (27) | size:6013 / node_N1114 (27) | size:6013.",
data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

levels(data_for_plot3$Species)= c("Culex pipiens"=make.italic("Culex pipiens"))

levels(data_for_plot3$Strain) <- c("Camping~Europe")

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

tiff("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_counts_serratia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_counts_serratia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_counts_serratia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_taxonomic_serratia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

img <- magick::image_read(paste0(path_serratia, "/HTML-OUTPUT/entropy.png"))
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

tiff("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_main_serratia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Dg_OLIGO_main_serratia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.1    openxlsx_4.2.3      Biostrings_2.54.0  
##  [4] XVector_0.26.0      IRanges_2.20.2      S4Vectors_0.24.4   
##  [7] BiocGenerics_0.32.0 cowplot_1.1.0       ggpubr_0.4.0       
## [10] RColorBrewer_1.1-2  forcats_0.5.0       stringr_1.4.0      
## [13] dplyr_1.0.2         purrr_0.3.4         readr_1.4.0        
## [16] tidyr_1.1.2         tibble_3.0.4        tidyverse_1.3.0    
## [19] ggplot2_3.3.2       phyloseq_1.30.0    
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0  ggsignif_0.6.0    ellipsis_0.3.1    rio_0.5.16       
##  [5] fs_1.5.0          rstudioapi_0.11   farver_2.0.3      fansi_0.4.1      
##  [9] lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  splines_3.6.3    
## [13] knitr_1.30        ade4_1.7-15       jsonlite_1.7.1    broom_0.7.2      
## [17] cluster_2.1.0     dbplyr_1.4.4      compiler_3.6.3    httr_1.4.2       
## [21] backports_1.1.10  assertthat_0.2.1  Matrix_1.2-18     cli_2.1.0        
## [25] htmltools_0.5.1.1 tools_3.6.3       igraph_1.2.6      gtable_0.3.0     
## [29] glue_1.4.2        reshape2_1.4.4    Rcpp_1.0.5        carData_3.0-4    
## [33] Biobase_2.46.0    cellranger_1.1.0  jquerylib_0.1.3   vctrs_0.3.4      
## [37] multtest_2.42.0   ape_5.4-1         nlme_3.1-149      iterators_1.0.13 
## [41] xfun_0.22         ps_1.4.0          rvest_0.3.6       lifecycle_0.2.0  
## [45] rstatix_0.6.0     zlibbioc_1.32.0   MASS_7.3-53       scales_1.1.1     
## [49] hms_0.5.3         biomformat_1.14.0 rhdf5_2.30.1      yaml_2.2.1       
## [53] curl_4.3          sass_0.3.1        stringi_1.5.3     highr_0.8        
## [57] foreach_1.5.1     permute_0.9-5     zip_2.1.1         rlang_0.4.10     
## [61] pkgconfig_2.0.3   evaluate_0.14     lattice_0.20-41   Rhdf5lib_1.8.0   
## [65] labeling_0.4.2    tidyselect_1.1.0  plyr_1.8.6        magrittr_1.5     
## [69] bookdown_0.22     R6_2.4.1          magick_2.5.2      generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       foreign_0.8-75   
## [77] withr_2.3.0       mgcv_1.8-33       survival_3.2-7    abind_1.4-5      
## [81] modelr_0.1.8      crayon_1.3.4      car_3.0-10        rmarkdown_2.7    
## [85] grid_3.6.3        readxl_1.3.1      data.table_1.13.2 blob_1.2.1       
## [89] vegan_2.5-6       rmdformats_1.0.2  reprex_0.3.0      digest_0.6.26    
## [93] webshot_0.5.2     munsell_0.5.0     viridisLite_0.3.0 bslib_0.2.4